// Copyright (C) 2011  Andrew H. Chan
//
// This program is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation, either version 3 of the License, or
// (at your option) any later version.
// 
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
// GNU General Public License for more details.

#include <iostream>
#include <iomanip>
#include <cstdlib>
#include <fstream>
#include <sstream>
#include <vector>
#include <stdexcept>
#include <assert.h>
#include<algorithm>

#include "LarssonSadakane.hpp"

#include "findMedian.hpp"

using namespace std;

typedef unsigned char uchar;
typedef unsigned int uint;
typedef unsigned long long uint64;

string LarssonSadakane::getSuffix(vector<uint64> const& T, uint64 index) {
    ostringstream iss;
    for(uint ii = index; ii<T.size(); ++ii) iss << T[ii] << ' ';
    iss << 0;
    return iss.str();
}

bool LarssonSadakane::suffixLessThan(vector<uint64> const &T, uint64 a, uint64 b) {
    if(a == b) return false;
    for(uint64 ii=0; a + ii < T.size() && b + ii < T.size(); ++ii)
        if(T[a+ii] < T[b+ii]) return true;
        else if(T[a+ii] > T[b+ii]) return false;
    if(a > b) return true;
    else return false;
}

void LarssonSadakane::qsort(uint64 h, vector<uint64> &V, vector<uint64> &SA, vector<bool> &L, uint64 index, uint64 length) {
    assert(length>0); // qsort should not be called with length 0 
    if(length < 7) { // selection sort if partition has fewer than 7 elements 
        uint64 pLow=index, pLast=index+length-1;
        while(pLow < pLast) {
            uint64 minKey = V[SA[pLow]+h];
            uint64 pEqual = pLow+1;
            for(uint64 ii=pLow+1; ii<=pLast; ++ii) { // search for minimum element
                uint64 tmpKey = V[SA[ii]+h];
                if(tmpKey < minKey) {
                    minKey = tmpKey;
                    swap(SA[ii], SA[pLow]); // swap to beginning of selection group 
                    pEqual=pLow+1; // new min found, so reset pEqual to point to one past end of min elements
                } else if(tmpKey == minKey) {
                    swap(SA[ii], SA[pEqual]); // swap to min elements group at beginning of selection group
                    ++pEqual;
                }
            }
            // update_group from pLow to pEqual-1
            uint64 groupNumber = pEqual-1;
            if(pEqual - pLow == 1) { // new sorted group
                V[SA[pLow]] = groupNumber;
                SA[pLow] = 1;
                L[pLow] = true;
            } else // new unsorted group
                for(uint64 ii=pLow; ii<pEqual; ++ii) V[SA[ii]] = groupNumber;
            // prepare to sort next selection group
            pLow = pEqual;
        }
        // finish updating group
        if(pLow == pLast) {
            V[SA[pLow]] = pLow;
            SA[pLow] = 1;
            L[pLow] = true;
        }
        return;
    } else { // quicksort if partition has 7 or greater elements
        // choose pivot
        uint64 pivot;
        if(length <= 50) { // if partition has 50 or fewer elements, select the median of first, middle and last elements
            uint64 a = V[SA[index]+h], b = V[SA[index+(length>>1)]+h], c = V[SA[index+length-1]+h];
            pivot = a > b ? (b > c ? b : (a > c ? c : a)) : (a > c ? a : (b > c ? c : b));
        } else { // if partition has more than 50 elements, select median of 5 evenly spaced elements
            uint64 incr = length>>2;
            pivot = findMedian5(V[SA[index]+h], V[SA[index+incr]+h], V[SA[index+incr+incr]+h], V[SA[index+incr+incr+incr]+h], V[SA[index+length-1]+h]);
        }
        uint64 pCurLow = index, pLow = index, pCurHigh = index+length-1, pHigh = index+length-1;
        uint64 curKey;
        while(true) {
            while(pCurLow <= pCurHigh && (curKey=V[SA[pCurLow]+h]) <= pivot) {
                if(curKey == pivot) swap(SA[pCurLow],SA[pLow++]); // swap equal elements to beginning
                ++pCurLow;
            }
            while(pCurHigh > pCurLow && (curKey=V[SA[pCurHigh]+h]) >= pivot) {
                if(curKey == pivot) swap(SA[pCurHigh],SA[pHigh--]); // swap equal elements to end
                --pCurHigh;
            }
            if(pCurLow>=pCurHigh) break;
            swap(SA[pCurLow++],SA[pCurHigh--]); // swap low with high
        }
        // pCurLow points one beyond the smaller elements
        // pCurHigh points one to the left of the larger elements
        pCurHigh = pCurLow-1;

        // swap equal elements from left toward center
        uint64 numEqLeft = pLow-index; // number of elements equal to pivot
        uint64 numLow = pCurLow-pLow; // number of elements less than pivot
        uint64 numEqRight = index+length-1-pHigh; // number of elements equal to pivot
        uint64 numHigh = pHigh-pCurHigh; // number of elements greater than pivot
        if(numEqLeft > numLow)
            for(uint64 ii=index, jj=pLow, numSwap=numLow; numSwap>0; --numSwap, ++ii, ++jj) swap(SA[ii], SA[jj]);
        else {
            for(uint64 ii=index, jj=pCurLow-numEqLeft, numSwap=numEqLeft; numSwap>0; --numSwap, ++ii, ++jj) swap(SA[ii], SA[jj]);
        }

        // swap equal elements from right toward center
        if(numEqRight > numHigh)
            for(uint64 ii=pCurLow, jj=index+length-numHigh, numSwap=numHigh; numSwap>0; --numSwap, ++ii, ++jj) swap(SA[ii], SA[jj]);
        else
            for(uint64 ii=pCurLow, jj=index+length-numEqRight, numSwap=numEqRight; numSwap>0; --numSwap, ++ii, ++jj) swap(SA[ii], SA[jj]);

        // sort smaller
        if(numLow>0) qsort(h, V, SA, L, index, numLow);

        // update group numbers
        // index is at beginning of smaller partition; pMiddle is at beginning of equal partition; pUpper is at beginning of larger partition
        uint64 pUpper = index+length-numHigh;
        uint64 pMiddle = index+numLow;
        uint64 groupNumber = pUpper-1; 
        if(pUpper - pMiddle == 1) { // new sorted group
            V[SA[pMiddle]] = groupNumber;
            SA[pMiddle] = 1;
            L[pMiddle] = true;
        } else // new unsorted group
            for(uint64 ii=pMiddle; ii<pUpper; ++ii) V[SA[ii]] = groupNumber;

        // sort larger
        if(numHigh>0) qsort(h, V, SA, L, pUpper, numHigh);
        return;
    }
}

int LarssonSadakane::verifyInverseSuffixArray(vector<uint64> const &T, vector<uint64> const &ISA, int printSuffixes) {
    // construct SA from ISA
    vector<uint64> SA(ISA.size());
    for(uint ii=0; ii<ISA.size(); ++ii) SA[ISA[ii]] = ii;
    if(printSuffixes) {
        cout << "Suffix array:" << endl;
        for(uint ii=0; ii<SA.size(); ++ii) cout << ii << ": " << getSuffix(T, SA[ii]) << endl;
    }
    cout << "Verifying suffix array..." << endl;
    bool badSuffixArray = false;
    for(uint ii=1; ii<SA.size(); ++ii) {
       if(!suffixLessThan(T, SA[ii-1], SA[ii])) {
           badSuffixArray = true;
       }
    }
    //if(badSuffixArray) throw logic_error("LarssonSadakane.cpp: Bad suffix array.");
    if(badSuffixArray) return 1;
    else return 0;
}

// V contain original sequence of uint64. V is overwritten with the inverse suffix array.
// Assumes V contains positive integers only, as 0 is used for the terminating character.
void LarssonSadakane::computeInverseSuffixArray(vector<uint64> &V) {
    V.resize(V.size()+1,0); // store intermediate keys; eventually holds ISA
    uint64 lenV = V.size(); // n+1
    vector<bool> L(lenV, false); // indicates if group is combined sorted group or unsorted group 
    vector<uint64> SA(lenV);// temporary storage 
    for(uint64 ii=0; ii<SA.size(); ++ii) SA[ii] = ii; //initialize SA
    
    // sort by first character
    qsort(0, V, SA, L, 0, SA.size());

    uint64 h = 1; // h-order
    while(!(SA[0]==lenV && L[0])) { // loop while not all elements are in one combined sorted group
        uint64 index = 0, lenSorted = 0;
        do {
            if(L[index]) { // if current index is a combined sorted group
                lenSorted += SA[index];
                index += SA[index];
            } else { // if current index is not a combined sorted group 
                if(lenSorted > 0) {
                    SA[index - lenSorted] = lenSorted; // combine sorted groups up to current index; store length in first index of combined sorted group
                    lenSorted = 0; // reset combined sorted group length
                }
                uint64 nextGroup = V[SA[index]] + 1;
                uint64 lenUnsorted = V[SA[index]] - index + 1; // length of unsorted group
                qsort(h, V, SA, L, index, lenUnsorted); // quicksort
                index = nextGroup; // move on to next group
            }
        } while(index < lenV);
        // if array ends with a sorted group, combine sorted groups up to end of SA
        if(lenSorted > 0) SA[index - lenSorted] = lenSorted;
        h *= 2; // double sorted prefix
    }
    // V contains ISA
    return;
}
